Preprocessing
Color Deconvolution
Separate stains
#Libraries
from matplotlib import pyplot as plt, patches
import numpy as np
#Enable nice output printing features
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'last_expr_or_assign'
import warnings
warnings.filterwarnings('ignore')
#Add other libraries as you see fit
import tifffile as tif
#Start coding here
#Import IHC image and split it to RGB
#IHC2
#img_ihc = tif.imread("IPQDA_23_ASS_E_DATA\IHC2.tif")
#H_E
img_ihc = tif.imread("IPQDA_23_ASS_E_DATA\H_E.tif")
#RGB channels
img_ihc_r = img_ihc[:, :, 0]
img_ihc_g = img_ihc[:, :, 1]
img_ihc_b = img_ihc[:, :, 2]
#Import cropped stain1, stain2 and background ROI images, OR import RGB vectors of the ROIs
#IHC2
#img_stain1 = tif.imread("IPQDA_23_ASS_E_DATA\IHC2_stain1.tif")
#img_stain2 = tif.imread("IPQDA_23_ASS_E_DATA\IHC2_stain2.tif")
#img_background = tif.imread("IPQDA_23_ASS_E_DATA\IHC2_back.tif")
#H_E
img_stain1 = tif.imread("IPQDA_23_ASS_E_DATA\H_E_stain1.tif")
img_stain2 = tif.imread("IPQDA_23_ASS_E_DATA\H_E_stain2.tif")
img_background = tif.imread("IPQDA_23_ASS_E_DATA\H_E_back.tif")
#End coding here
array([[[251, 241, 252],
[250, 240, 251],
[248, 238, 247],
...,
[255, 245, 255],
[255, 247, 255],
[255, 248, 255]],
[[252, 242, 253],
[247, 237, 248],
[248, 238, 247],
...,
[255, 246, 255],
[255, 247, 255],
[253, 243, 254]],
[[250, 240, 251],
[248, 238, 249],
[249, 239, 248],
...,
[255, 243, 253],
[253, 243, 254],
[250, 240, 251]],
...,
[[249, 237, 247],
[248, 236, 246],
[248, 236, 246],
...,
[251, 244, 251],
[251, 244, 252],
[252, 245, 253]],
[[253, 241, 253],
[252, 240, 250],
[250, 238, 248],
...,
[253, 246, 253],
[250, 243, 251],
[252, 245, 253]],
[[251, 240, 254],
[249, 239, 250],
[251, 241, 252],
...,
[255, 247, 255],
[252, 245, 252],
[254, 247, 254]]], dtype=uint8)
img_ihc.shape
(2309, 3877, 3)
#Inspect imported IHC image
plt.title("Raw IHC image")
plt.axis('off')
plt.imshow(img_ihc)
<matplotlib.image.AxesImage at 0x20280068550>
#Start coding here
#Calculate mean of image for each RGB channels. If you use RGB vectors, assign them directly to the variables here
mean_img_stain1 = np.mean(img_stain1, axis=(0, 1))
mean_img_stain2 = np.mean(img_stain2, axis=(0, 1))
mean_img_background = np.mean(img_background, axis=(0, 1))
#End coding here
print(mean_img_stain1)
print(mean_img_stain2)
print(mean_img_background)
[249.84285714 74.72857143 117.43571429] [66.8 13.70909091 71.57272727] [252.23859127 243.8735119 252.17361111]
#Convert RGB values to Hex color values for visualization
hex_img_stain1 = '#%02x%02x%02x' % tuple(mean_img_stain1.astype(int))
hex_img_stain2 = '#%02x%02x%02x' % tuple(mean_img_stain2.astype(int))
hex_img_background = '#%02x%02x%02x' % tuple(mean_img_background.astype(int))
print(hex_img_stain1)
print(hex_img_stain2)
print(hex_img_background)
#Visualization of RGB mean of cropped images
fig, axs = plt.subplots(1,3)
fig.suptitle('RGB mean of cropped images')
rectangle_stain1 = patches.Rectangle((0, 0), 1, 1, facecolor=hex_img_stain1)
rectangle_stain2 = patches.Rectangle((0, 0), 1, 1, facecolor=hex_img_stain2)
rectangle_background = patches.Rectangle((0, 0), 1, 1, facecolor=hex_img_background)
axs[0].add_patch(rectangle_stain1)
axs[1].add_patch(rectangle_stain2)
axs[2].add_patch(rectangle_background)
axs[0].set_title('Stain 1')
axs[1].set_title('Stain 2')
axs[2].set_title('Background')
axs[0].axis('off')
axs[1].axis('off')
axs[2].axis('off')
plt.show()
#f94a75 #420d47 #fcf3fc
#Calculate transmittances, T for each stain
T_stain1 = mean_img_background/mean_img_stain1
T_stain2 = mean_img_background/mean_img_stain2
OD_stain1 = np.log10(T_stain1)
OD_stain2 = np.log10(T_stain2)
print(OD_stain1)
print(OD_stain2)
[0.00414459 0.51367795 0.33189944] [0.57703507 1.25015598 0.54695207]
#Start coding here
#Normalize the absorbances
#Normalize by the squared of the sum of the vector values to the power of 2
OD_stain1_norm = OD_stain1 / np.max(OD_stain1)
OD_stain2_norm = OD_stain2 / np.max(OD_stain2)
#End coding here
print(OD_stain1_norm)
print(OD_stain2_norm)
[0.00806847 1. 0.64612359] [0.46157046 1. 0.43750706]
#Start coding here
#Combine OD_stain1_norm and OD_stain2_norm to form a normalized OD matrix M
M = np.column_stack((OD_stain1_norm, OD_stain2_norm))
#Calculate the deconvolution matrix according to Linear regression
MT = np.transpose(M)
MT_M = np.dot(MT, M)
inversed_MT_M = np.linalg.inv(MT_M)
D = np.dot(inversed_MT_M, MT)
D=D.T
#End coding here
print("M")
print(M)
print("M transposed")
print(MT)
print("Inversed M transposed multiplied with M")
print(inversed_MT_M)
print("Deconvolution matrix, D")
print(D)
M [[0.00806847 0.46157046] [1. 1. ] [0.64612359 0.43750706]] M transposed [[0.00806847 1. 0.64612359] [0.46157046 1. 0.43750706]] Inversed M transposed multiplied with M [[ 4.17951775 -3.82820822] [-3.82820822 4.21844559]] Deconvolution matrix, D [[-1.73326552 1.9162221 ] [ 0.35130953 0.39023737] [ 1.02561688 -0.6278959 ]]
#Convert pixel intensity to transmittance to absorbance according to Beer-Lambert Law on the IHC image
#Calculate the transmittance
T_img_ihc = mean_img_background/img_ihc
#Because of the logarithmic function in the next step, we assign all transmittance value less than 1 to 1
T_img_ihc[T_img_ihc<1] = 1
#Start coding here
#Calculate the absorbance
OD_img_ihc = np.log10(T_img_ihc)
#Coefficient matrix
coeffs = np.dot(OD_img_ihc, D)
#Extracting the individual coefficients from the coefficient matrix
#Which are essentially the orthogonal representation of the stains of the IHC image
coeff_stain1 = coeffs[:,:, 0]
coeff_stain2 = coeffs[:,:, 1]
#End coding here
print(coeff_stain1.shape)
print(coeff_stain2.shape)
(2309, 3877) (2309, 3877)
OD_img_ihc
array([[[2.88995293e-02, 4.52666181e-01, 2.74594839e-01],
[2.70631862e-02, 4.28123240e-01, 2.61820551e-01],
[1.79961663e-02, 3.91529438e-01, 2.40331635e-01],
...,
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[5.61218515e-03, 0.00000000e+00, 5.50029046e-03],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],
[[5.74192586e-02, 4.57745707e-01, 2.71365869e-01],
[3.81995524e-02, 3.91529438e-01, 2.31437922e-01],
[1.79961663e-02, 3.37946610e-01, 2.03042551e-01],
...,
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[2.13781077e-03, 0.00000000e+00, 2.99096775e-04],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],
[[5.93888514e-02, 4.37774626e-01, 2.58684837e-01],
[4.00836962e-02, 3.57780855e-01, 2.14178917e-01],
[2.16002905e-02, 3.11617671e-01, 1.86855790e-01],
...,
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],
...,
[[3.87152358e-03, 3.87164632e-01, 2.37346782e-01],
[3.87152358e-03, 3.70131293e-01, 2.22722690e-01],
[4.10991469e-04, 3.45771947e-01, 2.03042551e-01],
...,
[5.61218515e-03, 0.00000000e+00, 0.00000000e+00],
[1.08764251e-02, 0.00000000e+00, 3.75962888e-03],
[4.10991469e-04, 0.00000000e+00, 0.00000000e+00]],
[[5.61218515e-03, 3.78564461e-01, 2.34382303e-01],
[7.35985142e-03, 3.65975333e-01, 2.22722690e-01],
[5.61218515e-03, 3.49738134e-01, 2.08575039e-01],
...,
[3.87152358e-03, 0.00000000e+00, 0.00000000e+00],
[5.61218515e-03, 0.00000000e+00, 0.00000000e+00],
[3.87152358e-03, 0.00000000e+00, 0.00000000e+00]],
[[7.35985142e-03, 3.78564461e-01, 2.34382303e-01],
[7.35985142e-03, 3.61858767e-01, 2.19856050e-01],
[9.11457899e-03, 3.45771947e-01, 2.08575039e-01],
...,
[3.87152358e-03, 0.00000000e+00, 0.00000000e+00],
[2.13781077e-03, 0.00000000e+00, 0.00000000e+00],
[4.10991469e-04, 0.00000000e+00, 0.00000000e+00]]])
#Initialize the image absorbance container per stain
OD_img_ihc_stain1 = np.zeros((img_ihc.shape[0], img_ihc.shape[1], img_ihc.shape[2]))
OD_img_ihc_stain2 = np.zeros((img_ihc.shape[0], img_ihc.shape[1], img_ihc.shape[2]))
#Start coding here
print(OD_img_ihc_stain1.shape)
print(coeff_stain1.shape)
#Multiply the coefficients with the stain absorbance per stain. Do it independently for each RGB layer"""
"""for i in range(3):
OD_img_ihc_stain1[:, :, i] = OD_img_ihc[:, :, i] * coeff_stain1
OD_img_ihc_stain2[:, :, i] = OD_img_ihc[:, :, i] * coeff_stain2"""
OD_img_ihc_stain1 = np.reshape(coeff_stain1, img_ihc_r.shape)[:, :, np.newaxis] * OD_stain1
OD_img_ihc_stain2 = np.reshape(coeff_stain2, img_ihc_r.shape)[:, :, np.newaxis] * OD_stain2
#End coding here
(2309, 3877, 3) (2309, 3877)
array([[[0.03439603, 0.07451955, 0.03260283],
[0.03146729, 0.0681744 , 0.02982678],
[0.02098724, 0.0454692 , 0.0198931 ],
...,
[0. , 0. , 0. ],
[0.00421269, 0.00912687, 0.00399307],
[0. , 0. , 0. ]],
[[0.06824479, 0.14785347, 0.06468694],
[0.04654903, 0.10084924, 0.04412225],
[0.02243193, 0.04859914, 0.02126247],
...,
[0. , 0. , 0. ],
[0.00225547, 0.00488651, 0.00213788],
[0. , 0. , 0. ]],
[[0.07052009, 0.15278294, 0.06684362],
[0.0472861 , 0.10244611, 0.0448209 ],
[0.02635309, 0.05709441, 0.02497921],
...,
[0. , 0. , 0. ],
[0. , 0. , 0. ],
[0. , 0. , 0. ]],
...,
[[0.00546783, 0.01184615, 0.00518277],
[0.00693083, 0.01501575, 0.0065695 ],
[0.00474964, 0.01029016, 0.00450202],
...,
[0.00620555, 0.01344442, 0.00588203],
[0.01066418, 0.02310412, 0.01010822],
[0.00045444, 0.00098456, 0.00043075]],
[[0.00653002, 0.0141474 , 0.00618959],
[0.00985212, 0.02134479, 0.0093385 ],
[0.00938933, 0.02034213, 0.00889983],
...,
[0.00428085, 0.00927453, 0.00405767],
[0.00620555, 0.01344442, 0.00588203],
[0.00428085, 0.00927453, 0.00405767]],
[[0.00846246, 0.01833407, 0.00802128],
[0.00996379, 0.02158671, 0.00944434],
[0.01236891, 0.02679745, 0.01172407],
...,
[0.00428085, 0.00927453, 0.00405767],
[0.00236384, 0.00512129, 0.0022406 ],
[0.00045444, 0.00098456, 0.00043075]]])
#Convert absorbance to transmittance
T_img_ihc_stain1 = 10**(-OD_img_ihc_stain1)
T_img_ihc_stain2 = 10**(-OD_img_ihc_stain2)
array([[[0.92385534, 0.84232646, 0.92767781],
[0.93010657, 0.85472341, 0.9336266 ],
[0.95282416, 0.90059763, 0.95522769],
...,
[1. , 1. , 1. ],
[0.99034681, 0.97920389, 0.99084776],
[1. , 1. , 1. ]],
[[0.85458488, 0.71145351, 0.86161462],
[0.89836117, 0.79277649, 0.90339513],
[0.94965984, 0.89413039, 0.95222051],
...,
[1. , 1. , 1. ],
[0.99482006, 0.98881146, 0.99508944],
[1. , 1. , 1. ]],
[[0.85011936, 0.7034238 , 0.85734651],
[0.89683779, 0.78986685, 0.90194302],
[0.94112412, 0.87681019, 0.94410607],
...,
[1. , 1. , 1. ],
[1. , 1. , 1. ],
[1. , 1. , 1. ]],
...,
[[0.98748877, 0.97309188, 0.98813715],
[0.98416785, 0.96601585, 0.98498701],
[0.98912315, 0.97658452, 0.98968726],
...,
[0.9858128 , 0.96951734, 0.98654743],
[0.97574384, 0.94819111, 0.97699375],
[0.99895415, 0.99773553, 0.99900865]],
[[0.98507654, 0.96794928, 0.98584903],
[0.97757003, 0.95204004, 0.97872685],
[0.97861231, 0.95424056, 0.97971594],
...,
[0.9901914 , 0.97887102, 0.99070037],
[0.9858128 , 0.96951734, 0.98654743],
[0.9901914 , 0.97887102, 0.99070037]],
[[0.98070307, 0.95866292, 0.98169983],
[0.97731871, 0.95150985, 0.97848836],
[0.97192128, 0.94016169, 0.97336545],
...,
[0.9901914 , 0.97887102, 0.99070037],
[0.99457185, 0.98827705, 0.99485411],
[0.99895415, 0.99773553, 0.99900865]]])
#Clip each layer to 0,1
T_img_ihc_stain1[T_img_ihc_stain1>1] = 1
T_img_ihc_stain2[T_img_ihc_stain2>1] = 1
T_img_ihc_stain1[T_img_ihc_stain1<0] = 0
T_img_ihc_stain2[T_img_ihc_stain2<0] = 0
#Start coding here
T_img_ihc_stain1_norm = (T_img_ihc_stain1 * 255).astype(np.uint8)
T_img_ihc_stain2_norm = (T_img_ihc_stain2 * 255).astype(np.uint8)
#End coding here
array([[[235, 214, 236],
[237, 217, 238],
[242, 229, 243],
...,
[255, 255, 255],
[252, 249, 252],
[255, 255, 255]],
[[217, 181, 219],
[229, 202, 230],
[242, 228, 242],
...,
[255, 255, 255],
[253, 252, 253],
[255, 255, 255]],
[[216, 179, 218],
[228, 201, 229],
[239, 223, 240],
...,
[255, 255, 255],
[255, 255, 255],
[255, 255, 255]],
...,
[[251, 248, 251],
[250, 246, 251],
[252, 249, 252],
...,
[251, 247, 251],
[248, 241, 249],
[254, 254, 254]],
[[251, 246, 251],
[249, 242, 249],
[249, 243, 249],
...,
[252, 249, 252],
[251, 247, 251],
[252, 249, 252]],
[[250, 244, 250],
[249, 242, 249],
[247, 239, 248],
...,
[252, 249, 252],
[253, 252, 253],
[254, 254, 254]]], dtype=uint8)
#Display deconvolved image for stain 1
fig = plt.figure(dpi=300)
plt.title("Deconvolved image for stain 1")
plt.axis('off')
plt.imshow(T_img_ihc_stain1_norm)
fig.savefig('T_img_ihc_stain1_norm.tif')
#Display and export deconvolved image for stain 1
fig = plt.figure(dpi=300)
plt.title("Deconvolved image for stain 2")
plt.axis('off')
plt.imshow(T_img_ihc_stain2_norm)
fig.savefig('T_img_ihc_stain2_norm.tif')